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SUMMARY 


The efficiency and accuracy of several different techniques proposed for 
the numerical integration of ordinary differential equations with widely vary- 
ing damping properties are examined. The methods analyzed include the stan- 
dard Runge-Kutta method, a modified Runge-Kutta method proposed by Treanor, 
predictor -corrector schemes, and implicit procedures. First, the methods are 
studied as they apply to linear coupled differential equations with constant 
coefficients. The conclusions drawn from these studies are then tested in 
nonlinear cases by numerical calculations made for a gas in chemical non- 
equilibrium behind a normal shock wave. An implicit method is strongly 
recommended if the local eigenvalues are negative (damping) and widely 
separated. 


INTRODUCTION 


The numerical integration of the nonlinear equations arising in the 
study of a gas flowing in chemical nonequilibrium poses, in certain cases, a 
severe numerical stability problem. The existence of such a problem is well 
known (see refs. 1-5). It is also well known that numerical solutions involv- 
ing the use of standard Runge-Kutta and predictor -corrector methods (referred 
to as conventional methods) can be prohibitively expensive because they 
require excessive machine running time. 

Several attempts have been made to overcome the difficulty by Introducing 
numerical methods specifically designed to cope with the problem. These 
methods fall into two principal categories. In one, the nonlinear differ- 
ential equations are reduced to nonlinear difference equations, as in the con- 
ventional methods, but the coefficients in the differencing equations contain 
certain parameters which depend on the solution as it proceeds. If the 
differential equations were uncoupled, these parameters would be the local 
eigenvalues of the individual equations. A typical example of this class is 
the method proposed by Treanor. In the other category, the differential equa- 
tions are first locally linearized and the resulting linearized form is solved 
either exactly (ref. 6), approximately (ref. 7 ) ? or by finite difference 
methods (ref. 8). 



A basic objective which motivated the numerical research reported herein 
was to compare the efficiency of the various numerical methods available. To 
accomplish such an objective, it is essential that we know the fundamental 
relationship of a set of coupled differential equations and a corresponding 
set of coupled difference equations obtained from the former by the applica- 
tion of any given differencing scheme. To find such a relationship that 
applies rigorously to general nonlinear equations is quite outside the scope 
of this paper. However, the fundamental relations regarding the stability 
and accuracy of numerical solutions to coupled linear differential equations 
with constant coefficients can be rigorously defined. As a result of these 
definitions the numerical methods mentioned can be compared objectively as 
they apply to such linear equations. We then hypothesize that this comparison 
can be used to estimate their usefulness in nonlinear cases. 
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matrix of enclosed quantity 

inverse of matrix 

area 

matrix in locally linearized equations (see eq. (4)) 
determinant of enclosed quantity 
derivative of with respect to s 
value of F at step n 

effective distance a numerical method advances the integration after 
time for two evaluations of the derivatives 

enthalpy of the jth species 

unit matrix 

equilibrium constant for ith reaction 
number of equations 
number of chemical species 
Treanor parameter (see eq. (10)) 

production of species j in moles per unit volume per unit time 
universal gas constant 
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independent variable 

temperature 

velocity 

dependent variable In uncoupled equations (see eq. (5)) 
dependent variable in coupled equations (see eq* (2)) 
derivative of w with respect to s 
distance downstream of normal shock 

molar concentration of jth species, moles per unit mass 
upper bound of error permitted in numerical solution 
parasitic eigenvalue 
driving eigenvalue 
density 

degree of nonequilibrium of ith reaction 

Superscripts 

vector 

transpose of vector 
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THE BASIC NONLINEAR EQUATIONS 


The equations governing inviscid, one -dimensional, nonequilihrium flow- 
can "be written 
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Although the equations are written for general one -dimensional channel_flows, 
only the results for flow Behind a normal shock, for which A = 1 and A x = 0, 
are presented. This is done deliberately to isolate the basic numerical prob- 
lems discussed in the following sections. The fluid dynamic equations were 
differenced and numerically integrated just like the chemical equations even 
though they could have been integrated analytically. Thus, all the dependent 
variables are treated alike as they would be in a general two-dimensional flow. 


Equations (l) can be written in matrix -vector form by introducing the 
vector t?*, such that its_transpose w*^ = (u.p,T, 7 x , • • •> 7n) > ’ t ^ ie vector 
<?*, such that cf-*® = ( -puA x ,0,0,Q 1 , . . ., Q^) , and the matrix [B*]. The 
result is 


[B*] 


dw* 

dx 


c* 


In the general case the matrix [B*] can become singular (this occurs 
when u equals the sonic velocity) . Some of the numerical difficulties 
brought about by this occurrence can be avoided by introducing the new 
independent variable s such that 


x' 


& - <tet(B*) 
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and inverting the equation and multiplying "both sides by det(B*) to obtain 


v*' s = F* = det(B*)[B*]" 1 c* 
ds 


where det(B*) [B *] ~ ± is the adjoint of [B*]. Finally, we define the new 
vectors w and F with one element more than their starred counterparts 

= (w*Tjx) 


or 

w T = (u,p,T,7 1# . . .,7 n>x) 

and 


F T » [F* T ,det(B*)] 


This provides a set of m = N + 4 simultaneous equations 

£ - 


( 2 ) 


in which there is no explicit dependence of ? on the independent variable 
s. These are the basic nonlinear equations that motivated the numerical 
studies reported herein. 


LOCAL LINEARIZATION AND PARASITIC EIGENVALUES 


Consider the set of equations 


dw 

ds 


w' = F(tf) 


(3) 


These equations' are autonomous; that is, each equation is of the form 


w^ = F i (w 1 ,w 2 ,. . • ,w m ) 


where F^ has no explicit dependence on the independent variable s. If 
each F^ is expanded about a local point referenced as n* where s = nh, 


w i 


» F. + 
in 


(wi 



. + (w m - Wnn) 



+ 0[(w - w n ) 2 ] 
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Let the elements (a i( j) n of a matrix [A n ] be (dFq/dwj) n . Note that the 
neglected higher derivative terms involve (w - w n ) 2 . If w = w n+1 , this can 
be written h 2 [(w n+1 - w n )/h] 2 or h 2 (w^) 2 plus terms of 0(h 3 ) . Thus we find 

w' = [A n ]w + F n - [A n ]w n + 0(h 2 ) (4) 


Equation (4) gives, neglecting higher order terms, a locally linearized 
form of the original equations. Further, since the original equations were 
autonomous, the linearized equations have constant coefficients. It is well 
known that the complementary solution of equation (4) can be written 


^iji + C ij2 s + 




where K(j) is the multiplicity of the jth root to the characteristic equa- 
tion, J is the number of distinct roots, Aj, and m is the number of equa- 
tions. From a somewhat different point of view, the values of Aj (which can 
be complex) are equal in magnitude, sign, and multiplicity to the eigenvalues 
in the matrix [A n ] . 

Now it is typical of equations for nonequilibrium flow that in certain 
regions some of the eigenvalues are large, negative, real numbers and some 
are relatively much smaller in magnitude. Let these be represented by 
(A;]_) n = (^2) n = "I- 1 * respectively. When integrating equations with 

such eigenvalues, two cases can occur. One, we wish to resolve the effect of 
(Ai) n over the small region where e -1000 ^ Is significant. Then we must 
perform our calculations at points spaced very close together. Two, s is 
large enough for e -1000 ^ to be negligible compared to e”^ s . Then we need 
only use the much coarser spacing that resolves e -iJ,s . If the solution could 
be carried out analytically, this would pose no problem. However, if it is 
carried out numerically, using conventional Runge-Kutta or predictor -corrector 
schemes, violent instabilities occur for the coarse spacing. For this reason, 
when a set of differential equations has eigenvalues such as (Ai) n and (A2) n 
in the above example, we refer to those like (Ai) n as "parasitic" eigenvalues 
and those like (}\s) n as "driving" eigenvalues. Sets of equations having this 
property are often referred to as "stiff" equations. 

In order to discuss numerical methods used to solve equations with 
parasitic eigenvalues, we need the following brief discussion of some results 
from numerical analysis. 


A THEOREM REGARDING THE NUMERICAL SOLUTION 
OF ORDINARY DIFFERENTIAL EQUATIONS 

* 

Equation ( 4 ) represents a set of coupled, ordinary, differential 
equations with constant coefficients. In general, these can be uncoupled 
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insofar as they can he reduced to the Jordan canonical form (see ref. 9). For 
the purposes of numerical research it is often convenient to construct coupled 
equations from simple sets of uncoupled ones. Since this is the approach 
taken in appendix B, and since it illustrates the connection between coupled 
and uncoupled sets, we proceed to outline it here. 

Consider the special uncoupled set of linear equations with constant 
coefficients 


vi* = A1V1 + fi 

v 2 t = A2V2 + f 2 


These are special in the sense that their solutions have no terms of the form 
C'ijk^ with k > 0 in the complementary solution of equation (4 ) ; in other 
words, no coupled set formed from them will have a matrix with multiple 
eigenvalues. They can he written 


v* = [ L] v + f 


(5) 


where [ L] is a diagonal matrix (i.e., all off-diagonal elements are zero). 
Let 


w = [B]v 


( 6 ) 


where [B] is an arbitrary nonsingular matrix. Then, if 

[A] = [BHLKB]- 1 

and 

g = [B]? 


( 7 ) 


we can write 


v* = [A]w + g (8) 

where the eigenvalues of [L] and [A] are identical . Introduce the linear 
operator 2^ to represent any conventional 1 Runge-Kutta or predictor- 
corrector method with step size h. It is essential that, in a given step, 

1 More precisely, any set of linear difference -differential equations with 
constant coefficients applied in the same way to all of the simultaneous 
differential equations. All references to conventional numerical methods in 
this report assume a choice from this set. 
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the identical method be applied to all the equations. Then, if the numerical 
integration is started at the step n - j, at the step n we have 

v n = Z h (vj,Vj' ,f , j = n,n-l, . . . ,n-j 

The value of v n differs from the analytical value of v hy an amount equal 
to the error introduced by the numerical method. 

Consider the following theorem whose proof is given in reference 10. 

Theorem I : If v n = l h (vyv- ' ,£ j) and w n = Z h (wj ,gj) , then 

v n = [B]w n + er ; where er depends only on the round-off 
process used in the computation. 

The theorem can be paraphrased as follows: except for round-off error , solu- 

tions are the same whether or not a set of linear differential equations with 
constant coefficients is first uncoupled and then integrated numerically, or 
first integrated numerically and then uncoupled, no matter what conventional 
numerical method is employed, provided the same method is used on all of the 
equations. Two corollaries follow, both of which assume that machine limita- 
tions such as round off and floating overflow can be neglected. 

(1) Both the stability and accuracy of a conventional numerical method 
applied to equation (8) are independent of the magnitude and sign of the 
elements of [A] except as those elements determine the eigenvalues. 

(2) A step size that is chosen to resolve the effects of driving 
eigenvalues can be used to integrate equation (8) and the numbers obtained 
from the solution contain, to the desired accuracy, all information concerning 
the effects of the driving eigenvalues regardless of errors brought about by 
the parasitic ones. 

Corollary (2) is literally true even if the numerical method is unstable 
for the parasitic eigenvalues. However, the computing machine limitations 
neglected in formulating the corollary make this an impractical result for our 
purposes even if we were to uncouple the results. The important consequence 
of the corollaries is that to provide accurate and usable solutions to prob- 
lems with parasitic eigenvalues, we need only provide a method that is_ stable 
for large, negative, real values of Aih and accurate for any complex A 2 h 
relatively much smaller in magnitude. We next discuss some methods from this 
point of view. 


APPRAISAL OF METHODS 


With the preceding material as a background, we can now discuss the 
relative merits of several procedures used for the numerical integration of 
coupled equations with parasitic eigenvalues. These procedures include 
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explicit methods, implicit methods, and special methods that have been 
constructed for problems of this type. In order to compare such methods 
fairly, we first introduce a representative step size such that 

H = the effective distance a numerical method advances the integration 
after a time equal to the time required for two evaluations of the 
derivatives . 

Since in the numerical solution of nonequilibrium problems, the time required 
to calculate the derivatives is the predominant factor, H is the significant 
reference, rather than h, the step size used in the computation. Furthermore, 
we define the term real stability boundary as the largest real negative value 
of AH for which a numerical method is stable when used to integrate a set of 
coupled, linear, differential equations having A as an eigenvalue. 


Conventional Explicit Methods 


These methods encompass all standard predictor -corrector and Runge-Kutta 
techniques. About the simplest of such techniques, and one that has an 
accuracy that is acceptable for many practical cases, is given by 



v n+i 


v + hv T 
n n 


v, 


n 





(9) 


This method is referred to either as an Euler predictor with a modified 
Euler corrector or as a second-order, Runge-Kutta method. Actually, the 
method is unstable for problems with high-frequency, low-amplitude noise (see 
ref. 10), but this does not occur in the problems under consideration since 
the parasitic eigenvalues are real, not imaginary. Equations (9) have a 
truncation error led by (l/6)(AH) 3 and a real stability boundary equal to 
-2.0. Higher order Runge-Kutta methods are more accurate but less stable. 

For example, the third-order Runge-Kutta (Heune) method has a real stability 
boundary equal to about -1.6, and the standard fourth-order, Runge-Kutta 
method is limited still further to about -1.4. Most conventional predictor- 
corrector methods (Hamming's, Adams -Moulton, etc.) have real stability 
boundaries lying between -1 and 0. 


For all conventional explicit methods the step size is dictated by the 
largest eigenvalue of the system no matter whether It Is parasitic or driving. 


Treanor's Explicit Method 

Treanor's method (ref. 3) is explicit and was designed specifically to 
cope with the numerical integration of stiff equations. In constructing this 
method it is assumed that the basic equations can be approximated by a linear- 
ized form, but not the form given in equation (4). Instead, each equation is 
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( 10 ) 


differenced as if it could be written 2 

Vi’ = -(Pi) n Vi + (Ai) n + (Bi)nS + (Ci) n s 2 

where (P ) , (P 2 ) n , • • • are constants in a given step, but are allowed to 
vary from equation to equation, and Aj_, B^, and are all functions of 
only. The parameters (P i ) n would be the local eigenvalues of the individual 
equations if the equations were uncoupled. In practice they are formed by 
ratioing certain terms in the first two calculations made at the intermediate 
steps in a standard, fourth-order, Runge-Kutta process. 

Certain general statements can be made about this method. It reverts to 
the Runge-Kutta method if h is made small enough. Of course, if this prop- 
erty is used, the method will be no improvement on the conventional one. The 
method is excellent if the differential equations are "nearly" uncoupled. 
Unfortunately, however, it is not simple to make an a priori estimation of the 
degree of coupling in sets of equations. 

The most serious difficulty that can arise in the general use of 
Treanor 1 s method is discussed next. First, one can show (see appendix A) that 
if, in a given step, different differencing schemes are used on different 
equations in a coupled differential set, theorem I no longer holds. It fol- 
lows immediately that a method which permits Pj_ in equation (10) to vary 
with i will give solutions whose accuracy and stability will, in general, 
depend on the individual elements in [A n ] . Since, for a fixed set of eigen- 
values, these elements can have an extreme variation, such a dependency can 
lead to serious numerical errors . The actual values that the parameters 
(Pi) n are likely to take when found by Treanor’s method are studied in 
appendix B. It is shown that they are nearly always about the same and 
approximately equal in magnitude to the largest parasitic eigenvalue detect- 
able in w. However, the fact that this is not always the case is the 
principal deterrent to the general recommendation of the method. 

If the parameters in the differencing scheme proposed by Treanor 

are made to be the same for all i in a given step, theorem I does apply to 
the results. In such a case the method, when used to integrate coupled, 
linear, differential equations with constant coefficients, can be subjected 
to a complete analysis. This has been carried out and is presented in 
appendix C. The results regarding the stability for real negative eigenvalues 
are presented in figure 1. Since theorem I applies, the method acts on each 
eigenvalue as if the others were not present and a single figure gives the 
whole stability picture. If all the eigenvalues are such that -2 < A-h < 0, 
the Treanor differencing scheme with fixed P is stable for 0 < Ph < <». On 
the other hand, if an eigenvalue exists such that for it Ah < -2, a stability 
corridor starts to form, and the choice of P becomes critical. Two situa- 
tions can arise. First, if there is only one parasitic eigenvalue, or if all 


2 Another method which assumes this kind of local approximation is given 
in reference 11. However, because of the ambiguity in just how the values of 
(A i ) n , (B i ) n , and (Cq) n are to be formed in the case v’ = this 

method cannot be analyzed in detail. 
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Figure 1.- Stability "boundaries for Treanor's 
method for real negative A. 


the parasitic eigenvalues are 
clustered inside the corridor, the 
method is stable for indefinitely 
large separations between the para- 
sitic and driving eigenvalues, pro- 
vided P is appropriately chosen. 
Second, if two or more parasitic 
eigenvalues do not lie in the corri- 
dor, the method is still stable for 
coupled equations in which 
0 > Ajh > -10 if Ph is 8 . Since 
the method requires four evaluations 
of the derivatives in a single cal- 
culation step, this gives a real 
stability boundary AH equal to 
about 5 > an improvement of about 3*6 
over the standard, fourth -order, 
Runge-Kutta method. These observa- 
tions suggest certain modifications 
to Treanor's method if it is to be 
used in general; two such modifica- 
tions are discussed at the end of 
appendix B. 


Implicit Methods 

The methods discussed above are explicit. It has been known for some 
time that certain implicit numerical methods are unconditionally stable for 
all real negative values of Ah. Probably the most widely used of these is 
the Implicit modified Euler method described by 


n+i 


= v n + 


— h(v J 
2 s n+i 


+ v 1 


n' 


(ID 


which is even stable for all complex values of Ah with negative real parts. 
In studies of parabolic partial differential equations, this is commonly 
referred to as the Crank -Nicho Ison method (see, e.g., ref. 12, p. 2 6 k). 
Another method which also has second -order accuracy for small Ah, but is not 
so well known, is given by the two step equation 

v n+2 = 3 (Ki+i " v n + 2hv A+ 2 ) ( 12) 

Its use was suggested "by Curtiss and Hirschf elder (ref. 1, eq. ( 17 ) ) j and it 
is basically the method used in boundary -layer -theory studies by Davis and 
Flugge -Lotz (see ref. 13 ). Equation (12) is not quite as accurate as equa- 
tion (ll), but it is more stable, that is, the parasitic eigenvalues are more 
heavily damped. Methods that are unconditionally stable for real negative 
values of Ah and yet give higher order polynomial approximations for small 
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values of Ah are given in both references 1 and 8. They can he derived hy 
calculating the derivative v’ at n + k using the values of v at n + k, 
n + k - 1, . . . , n. The first five of these formulas are given in 
reference 14, pages 96 - 98 - 

The implicit method studied in this report is that given hy equation (ll). 
When applied to equation (4), there results the set of simultaneous equations 




(w n+1 - w n ) = hF n + 0(h 3 ) 


(13) 


where [A n ] and F n do not contain s explicitly, and F n = 

In order to compare this with other methods, it is necessary to estimate 
the time required to calculate the elements 8w^/dwj of [A n ] which, for 
nonlinear problems, must he reevaluated at each step. This time varies sig- 
nificantly with the problem and the details of the programming, hut a reason- 
able general estimate for problems with m simultaneous equations is to allow 
for m additional time intervals equal to that required for evaluating w^. 
Roughly, one more such interval is required to solve the simultaneous equa- 
tions, so the equivalent of about m + 2 calculations of w^ is required to 
advance the implicit method one step. Thus, approximately, h = H(m + 2)/2. 

The error of the implicit method is led by the term (Ah) 3 /l2 when applied to 
linear equations. For nonlinear equations the overall error is seen to 
remain at 0 ( h 3 ) . 

At this point the essential differences between the procedure proposed 
in this report and that proposed in reference 8 can be brought out. Basically, 
both processes depend upon an implicit method to insure stability. However, 
in reference 8 it is recommended that some of the equations (especially the 
fluid dynamic ones) be differenced by an explicit scheme. It was shown in the 
previous section that such a procedure invalidates the use of theorem I, which 
means that the stability and accuracy of the results become dependent on the 
size of the elements in the matrix. This can cause serious numerical diffi- 
culties. It is recommended here that all coupled equations be differenced at 
a given step by precisely the same differencing equation whether or not the 
method used is explicit or implicit. It is further suggested in reference 8 
that higher order implicit methods with step numbers greater than 1 can be 
used. In view of the error introduced into equation (13) by the nonlinear 
terms, the employment of methods embedding polynomials of order greater than 
2 cannot be justified in general. Further, the use of multistep methods 
detracts from the flexibility of the step-size adjustment available in one 
step methods. It is therefore proposed that when the use of an implicit 
method is warranted in the study of nonequilibrium flow problems, 
equation ( 13 ) is optimum. 


Locally Exact Methods 

Finally, let us consider the approach to the problem used in reference 6. 
In methods such as this the basic equations are linearized to the form 
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given by equation (4), and these linearized equations are then solved by 
actually finding the eigenvalues of the matrix [A n ] . Except for the fact that 
these eigenvalues must be determined numerically, this provides exact solu- 
tions to the linear equations, and we will refer to such methods as "locally 
exact." The essential differences in the detailed use of the various methods 
is summarized in the following table. For a fixed accuracy, proceeding 
downward in the table corresponds to increased computing times. 


Conventional predictor - 
corrector methods 

Conventional Runge -Kutta 3 4 
and Treanor methods 

Implicit methods 


Locally exact methods 


Must calculate Fj_ in equation (2) 
for i = 1, 2, 3> • • • , m 

Must, In addition, calculate 
dF^/chip for i = 1, 2, . . . , m 

Must, in addition, calculate all 
off-diagonal terms dF^/du^- for 
i,k = 1, 2, . . ., m; i / k 

Must, in addition, calculate all 
eigenvalues of [A n ] 


It is probably wise to estimate the largest local eigenvalue every once 
In a while in any event, but the exact calculation of all eigenvalues can be 
time-consuming, and, on this basis alone, 4 is to be avoided where possible. 

We seek, then, to find conditions under which the calculation of all the 
eigenvalues at every step can be avoided on the basis that such information 
increases neither the accuracy nor the stability of the results. In order to 
do this, consider the nature of our particular problem. Two possibilities 
occur. First, the linearized differential equations are themselves stable 
(inherent stability) . That is, the real parts of all eigenvalues in [A n ] are 
negative. Then, by theorem I, the implicit method given by equation (13) is 
just as stable as the locally exact methods, since the stability of both 
depends in the same way on exactly the same eigenvalues regardless of whether 
or not such eigenvalues are determined explicitly. Further, for general 
purposes, the implicit method is just as accurate as the locally exact methods 
since the order of error introduced by the linearization itself is as great as 
that caused by the differencing. Therefore, in the case of inherently stable 
differential equations, the continuous explicit calculation of the eigenvalues 
adds neither to the stability nor to the accuracy that is achieved by use of 
the implicit method given by equation ( 13 )* 

The other possibility mentioned above comes about when the linearized 
equations are themselves unstable. It Is conceivable in this case that the 

3 The standard, fourth -order, Runge -Kutta method is implied. For this 
and for Treanor* s method dF^/Sup can be determined by the calculations made 
at the Intermediate step. In the Runge -Kutta method it is not necessarily 
used although sometimes it is employed to monitor step size. 

4 The experience of the authors has been that standard subroutines that 
provide eigenvalues for general real matrices with parasitic eigenvalues are 
sometimes (which makes it worse) of questionable accuracy with regard to the 
driving ones. 
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locally exact methods would remain stable for significantly larger step sizes 
than would be possible for the implicit one. This condition can occur if the 
initial conditions are just those which require the coefficients of the 
unstable terms in the locally exact solution to be zero. The use of locally 
exact methods to solve nonlinear problems with this kind of inherent 
instability would require justification in individual cases. 


Recommendations 

On the basis of the above arguments, certain recommendations can be made 
which apply to the choice of a numerical method for integrating the equations 
(which are m in number, including those due to the fluid mechanics) govern- 
ing one -dimensional, nonequilibrium fluid flow. If is the largest para- 

sitic eigenvalue and | As | is the absolute value of the largest driving 
eigenvalue, and if, for accuracy, we require |A 2 |h<e, then, unless a special 
analysis of the particular problem indicates otherwise, it is strongly 
recommended that: 

(a) All m equations should be differenced by the same method for a 
given step, 

(b) The method should be the implicit one given by equation (ll) if 
-Ai/|A 2 | » (2/e)[(m + 2)/2] 

and two further , relatively weak 9 recommendations are: 

(c) The method can he the predictor -corrector combination given hy 
equation (9) if -Ai/|A 2 | < ( 2 /e)[(m -f 2 )/ 2] 

(d) The method can he the implicit one if > (2/e)[(m + 2)/2] 

Clearly it is not the purpose of this report to extol the implicit 
method as a universal technique to he used in solving nonequilihrium problems. 
No numerical method is better than all other methods in all respects, each 
having its own merits. In modern machine languages most methods of the kind 
considered herein are quite easy to code and represent only a small part of a 
sophisticated program. In general, it appears that the option of using 
more than one of them, as the occasion demands, appears to he the really 
optimum procedure. 


COMPARISONS OF NUMERICAL SOLUTIONS FOR THE FLOW OF AIR IN 
CHEMICAL NONEQUILIBRIUM BEHIND A NORMAL SHOCK WAVE 


All solutions were started hy applying the Rarikine-Hugoniot relations 
across a normal shock assuming a fixed chemical composition and local vibra- 
tional equilibrium for all species. The eight reactions shown in the insert 
in figure 2(h) were used. The rate equations and constants in them were the 
same as those used in reference 15 except that two reactions representing 0 
and N ionization are included. The recombination rates of the latter were 
assumed to he 100 times those given in reference 16. 
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Equations (l) were analyzed for a variety of cases } a representative one 
of which is shown in complete detail. For the representative case the free- 
stream conditions are 


T = 300° K 

P ot - 0. l6'95XlO -4 gm/cm 3 ) 
= 0.9144-X10 6 cm/sec 


(14) 


The degree of nonequilibrium, Xj_, for each reaction is defined as follows: 


*1 

1 - p7i 2 /r 4 Ki 

0 2 5 20 

x 2 

1 ~ P7 2 /V5K2 

N 2 5 2N 

X 3 

1 - P7 172 / 7 sK 3 

NO 5 N + 0 

x 4 

1 " 7176/7274^4 

N 4* "> NO 4* 0 

Xs 

1 - 727e/7i75K5 

0 4 N 2 ^ NO 4 N 

*6 

1 - P737s/7iK 6 

0 ^ 0 + + e" 

x 7 

1 - P 7379 / 72 K 7 

N % N + + e" 


1 - 737y/7 17 2^8 

F + 0 5 N0 + + e" 


where Kj_ is the equilibrium constant for the ith reaction and varies 

from unity initially to zero as equilibrium is approached. 


Numerical Solutions Using the Implicit Method 

The case described was computed over the range from just behind the shock 
to complete equilibrium, using the implicit method represented by equa- 
tion (13). (The first few steps were explicit in order to establish the 
matrix elements.) The results are shown in figure 2, the scale for which is 
in table I. It should be remarked that the elements in [A n ] were not found 
analytically. In every case they were calculated numerically by the equation 

SFj Fi(l.Olwj) ~ F i (0.99w -1 ) 

dwj ~ 0 . 02w j k 

which has an error 0 [(0.02wj) 2 ] . It should also be remarked that considerable 
computing time can be saved if the paths in the subroutine which calculates 
dFi/dwj are made different for u, p, T, and the species. As programmed, the 
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TABLE I . - VERTICAL SCALE PER MAJOR DIVISION IN 


FIGURES 2(a), 6, and 7 


Variable 

Scale 

cm/sec 

0.2X10 5 

p, gm/cm 3 

.4xio' 4 

p, dynes /cm 2 

.2X10 7 

T, °K 

. 4xio 4 

0 2 , moles/gm 

.1X10 -1 

N2, moles/gm 

.1X10 _1 

0, moles/gm 

.IXlO" 1 

N, moles/gm 

.IXlO -1 

NO, moles/gm 

.IX10 -1 

e”, moles/gm 

.1X10 -2 

N0 + , moles/gm 

.IXlO -2 

0 + , moles/gm 

.1X10 -2 

N + , moles/gm 

.1X10" 2 


machine computing time on an IBM 709^- averaged less than 90 seconds per case, 
and the total number of steps per case averaged less than 90. 

While use of the implicit method over most of the range is far from 
optimum on the basis of computing time, the results permit us to test the 
reliability of the local expansion procedures used in equations (4) and (15) 
when applied for many steps over a highly nonlinear region. This test can be 
made since the first part of the solution (over half of the steps) can be 
repeated using the same step intervals and conventional predictor -corrector 
methods which require no local expansions or construction of [A n ] . Such cal- 
culations also permit us to evaluate the local eigenvalues over the entire 
flow region so that the numerical results can be interpreted and correlated 
with the analysis presented in the previous sections. Finally, such results 
permit us to see how the integration of a set of 13, nonlinear, coupled equa- 
tions can proceed from a region having no parasitic eigenvalues into a region 
with several very large ones, without having stability difficulties. 

Before comparing the results of the implicit method with those calculated 
by explicit means, let us discuss the technique used to govern the step size 
as the calculations proceeded. 


The Method Used to Control Step Size 

Most techniques for adjusting the size of a step while the numerical 
integration is in progress depend upon the difference between a predicted and 
corrected value. One such method which tends to isolate the maximum negative 
eigenvalue in a set of stiff equations is presented in appendix B. When, 
implicit, methods are used, however, such a technique is not applicable since 




predicted and corrected values are not calculated. Furthermore, since the 
implicit methods of interest are stable for all negative eigenvalues, the max- 
imum eigenvalue is definitely not the one on which to base step size. As a 
result, the following method for controlling the step size of implicit methods 
was devised (it is not suitable for explicit methods when parasitic 
eigenvalues are present) : 

(1) After each step, compute Avf = ]w n+1 - w n | . 

(2) If any Awj < 0.001, ignore it in the following tests. (For the most 
part this limited the tests to the velocity and temperature variations.) 

( 3 ) If all Awj that passed test 2 were such that Awj/wj < 0.01, double 
the step size and take the next step. Otherwise, proceed to test 4. 

(4) If an^ Awj that passed test 2 was such that Aw^/wj > 0.10, halve 
the step size and take the next step. Otherwise, do not change the step size 
and proceed. 

No other tests (e.g., negative species concentrations) were made. This method 
gave the step-size variation shown by the symbols on the temperature curves in 
figure 2(a). 


The Eigenvalue Distributions 

A typical [A n ] matrix (step no. 59) is shown in table II. The diagonal 
elements are underlined and the eigenvalues (obtained from a subroutine listed 

in ref. 17) are also given. Because 



Figure 3*" Eigenvalues for solutions shown in 
figure 2. 


of limitations in the subroutine, the 
eigenvalues below 17.2 are not to be 
trusted (calculated eigenvalues were 
never used in any of the numerical 
integrations) . Such results were used 
to construct the curves shown in fig- 
ure 3, which gives the largest posi- 
tive eigenvalue and the largest 
negative eigenvalue throughout the 
region studied. The second and third 
largest negative eigenvalues over the 
latter part of the solutions are also 
included. The parasitic eigenvalues 
develop toward the end of the integra- 
tion. At first glance, the region 
over which the parasitic eigenvalues 
exist seems relatively small. But, 
since the plots are logarithmic, it is 
actually relatively very large. 


18 




TABLE II.- TYPICAL MATRIX [A n ] AND EIGENVALUES 



-0.716e 02 


0 

0.300E 01 
-0.318E-06 
-O.78OE-05 

0.411E-07 

0.802E-08 

0.379E-05 

O.236E-06 

0.719E-08 

O.583E-O7 

-0.241E-07 

0.986E-05 



-0.117E 11 
O.I87E 04 
0. 408E 05 
-0.908E 03 
-0.325E 02 
-0.197E 05 
-0.118E 04 
-0.374E 02 
-O.58OE 03 
-0.290E 03 
0.409E 04 


-0.162E-04 
0.266E-02 
0.172E-03 
0.214E-05 
-0.138E-02 
-0. 460E-04 
0.246E-04 
0.333E-04 
0.114E-03 
0.755E-04 



-0.860E 02 
0.607E 04 
-0.157E -00 
-0.487E 01 
0.20 IE 01 
0.344E 02 


Eigenvalues 


-0.123E 05 -0.782E 04 -0.204E 04 -0.948E 03 -0.124E 03 -O.806E 02 -0.576® 02 
0 0 0 0 0 0 0 


Matrix 


^Column 


-0.524E 09 
0.147E 01 
-0.125E 09 

-O.569E 02 
0.532E 03 
-0.139E 01 
-O.388E-OO 


0.599E 02 
-0.l6 IE -00 
-0.210E 01 
0.868E 00 
0.35IE 02 


9 

-0. 

10 9E 

10 

0. 

305E 

01 

-0. 

443E 

09 

0, 

16 5E 

04 

0. 

107E 

04 

6 

i 

243E 

02 

0 

266E 

03 

0, 

528E 

03 

-0 

.214E 

04 

-0 

•153E- 

-00 

-0 

. 412E 

02 

0 

■ 170E 

02 

0 

■ 335E 

02 


10 

0.134E 

11 

-0.375E 

02 

0 . 206E 

10 

0.783E 

04 

0.798E 

04 

-O.78OE 

04 

-0 . 422E ■ 

-00 

-0 . 102E 

03 

-0.173E 

02 

-0 . 778E 

04 

-0.415E 

02 

0.171E 

02 

0.411E 

02 


11 

12 

13 

0.124E 09 

0.119E 09 

-0 

-0.347E-00 

-0.332E-00 

0 

0.186E 08 

0.194E 08 

-0 

O.738E 02 

0.169E 02 

0 

0.162E 03 

0.204E 03 

0 

-0.571E 02 

-0.570E 02 

-0 

-O.358E-OO 

-0.343E-00 

-0 

-0-735E 02 

-O.665E 02 

-0 

-0.156E 02 

-0.151E 02 

-0 

-0.134E-00 

-0 . 121E-00 

-0 

-0.574E 02 

-0.100E 01 

-0 

0.458E-00 

-0.559E 02 

0 

0.293E 02 

0.264E 02 

0 


Eigenvalues 


0.10 IE 03 
0 


0.172E 02 0 

0 0 


-0.381E-04 

0 
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Figure 4 is a plot of Ah in which A is the largest positive 
eigenvalue. The product | Ah | max in which A is the largest negative eigen- 
value is shown in figure 5* These clearly establish the region outside of 
which numerical methods with limited stability boundaries would fail. For 
example, the second -order Runge-Kutta method (eq. (9)) would be unstable when 
|Ah| m ax exceeded 2. 



icr 6 icr 5 icr 4 io ~ 3 i0' z i0"' 


x, cm 

Figure 4.- Product of h and largest positive eigenvalue. 

Results Using Other Methods 

The initial values given in equations (l4) provided a model on which to 
base further numerical tests using different methods. First, the explicit 
second -order Runge-Kutta method given by equation (9) was combined with the 
method for step control given above and used to integrate the basic formulas 
directly in the form presented in equation (3)* In the range where |Ah| max 
remained less than 2, the results were almost identical to those given by the 
implicit method applied to equation (4). This occurred for values of x less 
than about 0.002 (see fig. 5)- A comparison of the two results is shown in 
table III. In each case, about 62 steps were taken and the variation of step 
size was essentially the same. (The slight difference in s is accounted for 
by the few, fixed, explicit steps initially used to develop a reliable [A n ] 
matrix before the implicit method was turned on.) The results show that 
numerical methods that make use of 

(a) Local expansion of nonlinear coupled equations, 

(b) Numerical evaluation of the elements in [A n ] , and 
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x, cm 

Figure 5*- Product of h and largest negative eigenvalue. 


TABLE III . - COMPARISON BETWEEN RESULTS OF NUMERICAL INTEGRATIONS 
AFTER ABOUT 62 STEPS USING EXPLICIT AND IMPLICIT METHODS 



Explicit 
(eq. (9)) 

Implicit 
(eq. (11)) 

s 

0.20140E-02 

0.20100E-02 

X 

.18588E-02 

.18560E-02 

u 

.73188E 05 

•73179E 05 

p 

.21177E-03 

.2118 IE -03 

T 

.12326E 05 

.12327E 05 

0 

.13125E-01 

.13122E-01 

N 

.3018 IE -01 

.30177E-01 

e " 

.25602E-02 

.25616E-02 

0 2 

.10907E-04 

.10903E-04 

n 2 

.11215E-01 

.11218E-01 

NO 

.4879LE-03 

.48 814E-03 

NO + 

.8I259E-05 

.81074E-05 

0 + 

.10 16 IE -02 

.10196E-02 

N+ 

. 15359E-02 

•15339E-02 
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(c) Implicit integration of the equations thus derived 


give an excellent numerical solution as far as accuracy is concerned, even 
when the spread of elements in [A n ] is as drastic as that shown in table II. 

It should be remarked that all calculations were made in 8 -place, 
floating-point arithmetic. 

When the explicit method was used in the numerical integration for 
x < 0.002 and the implicit method for x > 0.002 (at all times monitoring the 
step size by the method for step control given above), the results at x = 3*24 
were even closer to the "pure" implicit method than the comparison shown in 
table III. A plot of the complete results of these computations is shown in 
figure 6. The total machine time required for the combined explicit -implicit 
calculation was less than half of that required to calculate the entire range 
using the implicit method. 



x, cm 


Figure 6.- Numerical calculation using the explicit method for x < 0.002 (62 steps) and 
implicit method for x > 0.002 (29 steps). Step control given on page 18. 

Another calculation was made using the explicit method given by 
equation ( 9 ) for x < 0.002, and then switching to Treanor’s method for 
x > 0.002. The results are shown by the solid lines in figure 7- The method 
for step control given above was used throughout and Pq was set to zero if 
it was calculated to be negative. The computation was continued for the same 
number of steps ( 29 ) as that used for the implicit continuation shown in fig- 
ure 6. Under the imposed conditions, the results are what would be expected 
for the parasitic eigenvalue spread shown in figure 3 and the step control 
employed. 
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A final calculation, made using the computer program described in 
reference 15, is shown by the dashed curves in figure 7* For the lower values 
of x these lines coincide with the solid lines. The computer program of 
reference 15 makes use of Treanor* s method but with a suitably designed step 
control described in the reference. The program was made to compute for 
5 minutes and progressed to a value of x equal to O.W>, at which point the 
step size was 0 . 002 . At every computed point the results agreed with those 
calculated by the implicit method and there were no visable oscillations. 

This shows that, with proper step control, Treanor 1 s method can give excellent 
results. Its limitation in carrying solutions into regions with large para- 
sitic eigenvalues is demonstrated by the fact that its running time was six 
times that required for the calculations shown in figure 6. 



10" 5 I0" 4 JO" 3 I0“ 2 10“' 


x, cm 

(a) Thermodynamic variables. 

Figure 7*- Numerical calculation using equation (9) for x < 0.002 (62 steps) and Treanor ? s 
method - with the step control given on page 18 - for x > 0.002 (29 steps), solid lines. 
Numerical calculation using Treanor's method as programmed in reference 15> dashed lines. 
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(b) Molecules and atoms. 



x, cm 


(c) Electrons and ions. 
Figure 7 *“ Concluded. 
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CONCLUDING REMARKS 


If a set of ordinary differential equations can be linearized with regard 
to the dependent variables for a given value of the independent variable, then 
local behavior can be related to the eigenvalues of the matrix constructed 
from the linearized form. Then, if the differential equations are integrated 
numerically in a proper way, the stability (and accuracy, if the equations are 
autonomous) of the integration depends only upon these same eigenvalues and 
is independent of the detailed coupling of the set. When some of these eigen- 
values are parasitic, conventional explicit methods can require excessive 
machine computing times. If properly applied, specially designed explicit 
methods, such as that due to Treanor, can give considerable reductions in com- 
puting time. However, we conclude from the results of this study, that if 
more than one very large parasitic eigenvalue occurs, the differential equa- 
tions should be linearized locally and integrated implicitly (using the 
simplest implicit method compatible with the desired accuracy). 


Ames Research Center 

National Aeronautics and Space Administration 

Moffett Field, Calif., 94035.* May 31, 1967 
129 - 04 - 03 - 02 - 00-21 
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APPENDIX A 


EFFECT OF USING SEVERAL DIFFERENCING SCHEMES IN ONE STEP 


I n this appendix it is shown that when, for a given step, different 
schemes are used to difference the various differential equations in a coupled 
set, the stability and accuracy of the result will generally depend upon the 
individual elements in the [A] matrix of the differential set. 


Consider the two linear differential equations 

w’ = anw + a 12 v 
v* = a 2 iw + a 22 v 

Apply the Euler predictor 


u n+i = u n + 


to the first, and the implicit modified Euler method 


u 


n+i 


u n + p k( u n+i + u n'^ 


(Al) 


(A2) 


(A3) 


to the second. There results 

w n+i = w n + h ( a n w n + a i 2 ^ n ) 

v n+i = v n + 2 h ( &2lW n+i + a 22 v n+i + a 2i ¥ n + a 22 V n ) 

Introduce the operator E = exp[h(d/dx)] and combine terms. We find the 
matrix equation 


E - (l + a u h) 
| a 21 h(E + 1) 


E 


-a 12 h 


V 

n 

1 ~ 2 a 22 (F + l)_ 


— n _ 


having a characteristic equation which can be written 

4li - S x (E,h) a 12 

a 2 l 3*22 “^2 (B,h) 

in which S x (E,h) and S 2 (E,h) are the operators 


det 


= 0 


(At) 


2 6 



E - 1 
h 


Si = 


S 2 = 


2(E - 1) 
h(E + 1) 


(A5) 


It is at once evident that the form of results entirely from the choice 

of an Euler predictor for the first row, and the form of S 2 is "brought about 
entirely by the use of an implicit modified Euler method for the second row. 
The generalization is clear. If a variety of methods are used on the indi- 
vidual equations, a variety of operators appear on the diagonal of [A]. 

Returning to the simple example given by equation (A^) , and using the 
identities 


Al + ?\2 = a ll + a 12 


A 1 A 2 - a lia*22 ~ a 12 a 21 


where Ai and A 2 are the eigenvalues of 

a ll a 12 

a 21 a 22 

we can reduce equation (Ak) to 

S1S2 - (Ai + ^2)^2 + a 22 ( S2 * Si) + A1A2 = 0 


(A6) 


Suppose the differencing equations (A2) and (A3) had been the same. Then 
Si = S 2 = S and the left side of equation (A 6) factors into (S - Ai)(S - A 2 ) 
which shows no dependence on any of the elements, a-jj, regardless of the 
functional dependence of S on E and h. Obviously this would be true in 
general. However, for the different values of S L and S 2 in the example, 
equation (A 6) becomes 

h 2 

+ — AiA 2 + ha 22 


E" 


ha 


1 - 


22 


+ E 


-2 - h(Ai + Aa) 


1 + h(Ai+ A 2 )+ A 1 A 2 “ 


ha 


*22 


= 0 (A7) 


which clearly depends on the element a 22 as well as on the eigenvalues. 

As a simple illustration, let Ai = 1, A 2 = -1, and t = (l/2)ha 22 . The 
roots to equation (A7) are then 


E = 1 ± h 



+ 0(h 2 ) 
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The two matrices 







have the same eigenvalues (+1 and -l) , but have values of T equal in one 
case to -8 • 5h and to +8. 5b in the other. Use of the mixed differencing 
schemes in the example is clearly inadvisable. 
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APPENDIX B 


ANALYSIS OF THE PARAMETER P IN TREANOR’ S METHOD 


The ratio, P., between the difference in the predicted and corrected 
values of the derivative and the function calculated at the intermediate step 
in a standard, fourth-order, Runge-Kutta method is a basic parameter used in 
the method proposed by Treanor. In this appendix we analyze the significance 
of Pq when it is calculated from linear equations with constant coefficients. 
A set of "representative" stiff equations is then introduced and analyzed by a 
variety of methods. It is shown that for any method an effective Pq (that is, 
a number equal to that which would have resulted if the Runge-Kutta method had 
been used) can be calculated as the integration proceeds. The consequence of 
the results with regard to the accuracy and possible modification of Treanor 1 s 
method is discussed and a variety of examples are shown. 


The predictor -corrector equivalent of the standard, fourth -order, 
Runge-Kutta method is 


(1) 

V n+(i/2) " V n + 


i\ 
2 ) 


v„ 




( 2 ) /h\ (i)' 

V(i/2) = V n + ( 2 ) V(i 72 ) 


• (3) . V + 2 ^ V (2) ' 


n+i 


n 


2 ) n+(i/2) 


> (Bl) 


1 f h 

v n+i v n + 3 l 2 


A 3 )' + 21 V 
n+i 


(2)< (l)» 

n+(i/2) + V(i/2) 


+ V, 


n 


where the superscript designates a family generated in a single step. Let v 
represent the ith element in the vector w. In Treanor’ s method the value 
of Pq used in equation (lO) is determined by the equation 


(2)’ (i)» 

V(i/ g ) - V n + (!/2) 
"Ta) (1) 

n+(i/2) n+(i/2) 


(B2) 


Consider the coupled set of linear differential equations with constant 
coefficients 




dw 

ds 


= [A]w( s) + f(s) 


(B3) 
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which are not autonomous, and in which [A] is a matrix of constants. The 
explicit results of applying equations (Bl) to (B3) are, for the first family. 


"L(l/2) (j 1 ' + (a) ^*-0 + (2 ) 


- ([A] + (If [A] 2- ) \ + (|) [A]f n + f n+1 


n+( 1 / 2 ) 


(B4) 


and, for the second family, 


w 


( 2 ) 


n+ 


h 


(1/2) + (2) ^ + (2) ^ ' \ 2 J JJ ‘ n ' \ 2 ) "n+i 


+ (£) [A]f n + (SU 


(a)« 

W n+(x/2) 


= ( U] + (£) [A] 2 + (^ [A] 3 ^ n +( r |f[A]^r n + (|) [A]f n+1 4- f. 


n \ 2 


n+i n+i 


(B5) 


From these we form the ratio 


-> 

P = 
n 


[A] j[A] 2 w n + [A]f n + ' n+ ^ ' /^ " ' ^ 
{[A]^ n + [A]f n + 


(B6) 


where the division is defined to mean that an element in the vector in the 
numerator is divided by the corresponding element of the vector in the 
denominator. 

In studying the effect of parasitic eigenvalues, as we have defined 
them, we can consider the dominant eigenvalues of [A] to be real and distinct. 
Now let g be any linear combination of the eigenvectors of [A]. A well- 
known result under these conditions (see, e.g., ref. 9 , pp. 205 and 206) is 
that if 

A = (B 7 ) 

[A] 3 & 

— >• 

and division is defined as above, when j is increased, all elements in A 
approach the dominant eigenvalue whose eigenvector is given a nonzero weight 
in the formulation of g. 
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We now seek to find how closely the results just presented apply to the 
value of ? n as n is advanced in a numerical solution. Our procedure is to 
make a simple numerical experiment starting with three uncoupled equations and 
their initial values such that 


Vi* = -Vi + 1 , 

vi(0) = 0 

v^ 1 = ” 500 v ,2 y 

v 2 (0) = 0.002 

V 3 1 = - 1000 V 3 + 1 , 

v 3 (0) = 0.001 


(B8) 



S 


The exact solutions are shown in 
sketch (a) , and the value of [ L ] as 
used in equation (5) is given by 


[L] = [I] 


-1 

-500 

-1000 



Sketch (a).- Exact solutions of equations (B8). 


The eigenvalue -1000 is everywhere 
parasitic since v 3 does not vary at 
all, and the eigenvalue -500 quickly 
becomes parasitic relative to -1. 

These equations are coupled by a set 
of transformations given in equa- 
tions (6) and (7) the matrix used for 
[ B ] (chosen at random) is 


[B] = 


'- 1.8944 

-.14343 

1.4721 


-0 .84661 

- . 1^809 
- 1.1597 


- 1.2519 " 
-.40617 
-.41521 


The corresponding matrix [A], which equals [B][L][B] is 1 


[A] 


-19409 
-46 57-6 
314.68 


74820 -17686 

17717 -4267.3 

-2187.2 190.85 


(B9) 


(B 10 ) 


The coupled equations for w take the form 


X A11 computations were made in double precision to isolate roundoff 
errors. We present [A] only to illustrate the spread in the values of the 
elements, and both [A] and [B] have been truncated, for presentation, to five 
places . 
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v» = [A]w + [B]f 


(Bll) 


The behavior of the exact solutions 
for w is shown in sketch (b) . 

Several numerical methods were used to 
integrate equation (Bll) for about 
100 steps , and^for each method the 
variation of P was recorded. 


Before discussing the results, it 
is necessary to say something about the 
step-size control used in the integra- 
tion. In cases with parasitic eigen- 
values, the magnitude of the largest 
negative eigenvalue is not necessarily 
the best parameter with which to judge 
step size. A variety of methods for 
controlling step size has been pro- 
posed, but we do not wish to list or 
judge them here. The particular procedure -adopted in all the numerical 
methods discussed in this appendix is the following: 

(1) Choose an initial step size. ^ 

(2) Advance one step and compute 

Aw = max|w n+1 - w n | 

( 3 ) If Aw > 0.05 halve the step size and repeat. 

(t) If Aw < 0.0125 double the step size and repeat unless 

(a) An imposed maximum has been reached, 

(b) The step size was previously doubled in this same step, j 

The study of equation (Bll) is a controlled experiment and this procedure was 
designed especially for it. The procedure is not recommended in general. 

Using the implicit method given by equation (ll) together with the step- 
size control just described, the equations for w were integrated for 100 
steps starting with a step size equal to 0.02. The results, plotted in 
sketch (b), were excellent for the full 100 steps. The value of s was 
advanced to 7.74- (the symbols show the actual step locations) and the final 
Ik step sizes were equal to 0.32, the maximum allowed. When, after 100 steps, 
w was uncoupled (using the relation v = [B]" 1 ^), it compared to the exact 
analytical solution for v as follows: 


> (B12) 



Sketch (b).- Variation of v in equation (Bll) 
calculated by the implicit method using step 
control (B12). 



Exact 

Numerical 

Vl 

v 2 

V 3 

0.999565 
. 000000 
.001000 

0.999582 

.000000 

.001000 
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Under the constraints imposed hy the method used for step control, the 
step size remained at 0.02 for the first 55 steps. During the first 10 of 
these, P was essentially the same for all three w T s and equal to about 500. 
A transition occurred between steps 10 and 40 during which wild fluctuations 2 
of P were noted. From step bO through step 73 all values of P were again 
very uniform and all remained equal to about 1. At this point the step size 
was doubled from 0.04 to 0.08 and, as further doubling occurred and the w T s 
approached their asymptotic values, ? began to fluctuate. The actual 
variation of P is shown in sketch (c) . 
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Sketch (c).- Value of P3 calculated from 


With the result given by (B6) 
and the discussion concerning (B7), 
the variation of ? shown in 
sketch (c) can be explained. In the 
first place we notice that the 
dominant eigenvalue (the one associ- 
ated with V3 and equal to -1000) 
was never detected by equation (B 6) . 
The reason is simple. The exact solu- 
tion for V3 under the given initial 
conditions is a constant. In the 
terminology used to discuss equa- 
tion (B7) , the corresponding eigen- 
vector has a zero weight in the 
construction of w. To be sure trun- 
cation errors will begin to excite 
v 3 but, since the method is abso- 


equation (b6) using results from implicit 
method. Variations of Pi and P2 were 
nearly identical with P3. 

lutely stable for all negative A, this excitation is continually suppressed 
and has no effect on the value of P throughout most of the calculation. On 
the other hand, the eigenvector associated with v 2 is activated by the ini- 
tial conditions and, for a while, appears in all three elements of P. But 
relative to Vi, V2 is very heavily damped and in any numerical integration 
that is stable for all negative A, its influence eventually disappears from 
the leading significant figures that affect P. For the implicit method used, 
this disappearance is completed when s ^ 0.6 and, at that point, the largest 
eigenvalue detectable in w is -1. The disappearance of v 2 does not occur 
at once, though, and for many steps a "struggle for dominance" (ref. 9, 
p. 206) between the eigenvalues associated with v x and v 2 takes place. This 
phenomenon is quite evident in sketch (c) for values of s around 0.3* Even- 
tually all elements of w approach a constant and, although the solution 
itself is quite stable, P begins a meaningless fluctuation brought about by 
small truncation errors in all of the terms. 


Next, equation (Bll) was analyzed by a standard, fourth-order Runge- 
Kutta method using the identical step -size control. The results for are 
shown in sketch (d). The method is unstable for | Ah| > 2.8 (|ah| > 1.4) so, 


o —r 

It should be emphasized that P was in no way used in the numerical 
integration. 
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Sketch (d).- Solution of equations (Bll) using 
standard, fourth-order, Runge-Kutta method 
with step-size control (B12) . 


!500 r 


with this kind of step control, it 
automatically settled 3 on a step size 
equal to 0.0025 for most of the inte- 
gration. The scallop effect (which 
can also he noticed in the Runge-Kutta 
results presented in ref. 3) occurs 
because every once in a while , for 
just one step, use of the unstable 
step size, 0.005, is permitted by the 
particular test employed. This causes 
a relatively large error in due to 
the presence of V3. Immediately the 
method is forced to return to its 
stable step size and remain there 
until this error is sufficiently 
reduced, at which point the phenomenon 
is repeated. The value of s is 
advanced to 0.315 in 100 steps, and 
for this s the uncoupled values 4 of 
v are 
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Sketch (e).- Value of P 3 calculated from 

equation (b 6 ) or equation (B2) using results 
of Runge-Kutta method. Variations of Pi 
and P 2 are nearly identical with P 3 . 


The value of P generated by 
this method (once again, this value 
was never used in the integration) is 
shown in sketch (e) . The method 
started itself with step sizes varying 
between 0.02 and 0.005 because it 
could not immediately detect the 
largest negative eigenvalue. After 
13 steps, however, the error intro- 
duced by the unstable integration 
brought the presence of V3 clearly 
into view, and throughout the res^b of 
the calculation all elements of P 
were about 1000. 


Next, Treanor f s method with the same step-size control was used to 
integrate equation (Bll) . This time, of course, P was both calculated and 
used at each step. The results for ^ are shown in sketch (f). After the 
third step the method had IT found M the largest eigenvalue and all three 

3 The step-size control being used is quite inefficient for this and other 
methods with limited stability boundaries. In the presence of parasitic 
eigenvalues such methods increase their step size until they reach their 
stability boundaries and then proceed to first double and then halve at almost 
every step. 

4 If the uncoupling is made just after an unstable step size is used, the 
error in v 3 is about five times that shown. 




Sketch (f ) . - Solution of equations (Bll) 
using Treanor's method with step-size 
control (B12) . 


elements of P were nearly 1000 for 
steps 3 through 5. In this range a 
step size of 0.02 passed all tests. 
This procedure is unstable for the 
-500 eigenvalue, however, and at the 
sixth step Pi suddenly became nega- 
tive, while P 2 and P3 changed to 97 
and 364, respectively. The rule in 
using Treanor’s method is to set nega- 
tive values of Pj equal to zero and 
proceed. This can cause an abrupt 
change in the accuracy and stability 
of the integration, a change which, 
in general, depends critically ( see 
appendix A) on the magnitude and sign 
of the individual elements in [A]. In 
the present example the result was the 
same as If a completely new set of 
initial values were introduced, and 
the Integration departed drastically 
from the desired solution. At each of 
the three abrupt vertical rises in w 3 
shown in sketch (f), the sign of one 
of the elements in ? suddenly became 
negative while the other two remained 
at large positive values. It was fur- 
ther noted that a few steps before 
each of the less abrupt, nearly hori- 
zontal turns the same behavior in P 
occurred. 


To give Treanor’s method a chance, it was again used to integrate 
equation (Bll), but with a step size fixed at 0.008, and all. elements of P 
fixed at 1000. According to the stability boundaries shown in. figure 1, this 
integration should be stable throughout. After 100 steps s had been 
advanced to 0.8 and the uncoupled values of v were 



Exact 

Numerical 

Vi 

v 3 

0.550671 

.000000 

.001000 

0.550576 

.000001 

.001000 


To verify further the stability boundaries in figure 1, this. type of run was 
repeated with all elements of P again fixed at 1000, but. with a fixed step 
size equal to 0.01. The results were an interesting verification of theorem 1 
and corollary 2. After 50 steps s was advanced to 0.5 and when w was 
uncoupled, there resulted for v 
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.000000 
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All elements of the computed 5 P (which, of course, was not used) were almost 
exactly 500 throughout the run. Notice that, although an enormous error is 
made in v 2 ( all elements of w were greater in magnitude than 10 s at the 
fiftieth step), the two uncoupled terms for which the method is stable 
remained accurate through the fourth significant figure . 

The results of the test runs made using Treanor's method with fixed 
values of P. and step size suggested the following modifications to his 
method whicn would make it more reliable in general cases. One modification 
would be to 

(1) Use the same differencing scheme 

(2) Find P by means of equation (B2) and let P max = max(P) 

(3) If P max is negative, set P max = 0 

(4) Reset all elements of P with Pmax so that all equations are 
differenced by the same equations in a given step. 

When applied to equations with widely different parasitic eigenvalues, this 
method has a real upper stability boundary 1 7\h] = 10 and |AH| = 5* 

This modified method was combined with the same step control employed in 
the previous examples and used to integrate equation (Bll) . The results for 
w are shown in sketch (g) . Since it was limited by the unrefined doubling- 
halving step control procedure, the method was forced to fluctuate between 
step sizes of 0.005 and 0.01 for which it was stable and unstable, respec- 
tively. The scalloped effect caused 
when this kind of step control is used 
in methods with bounded stability is 
evident. The computed values of ? 
are shown in sketch (h) , In each step 
all three values of P were made to 
equal the maximum one shown. The 
method advanced s to 0.935 in 100 
steps. Define the error terms 


Av j “ exact " ( v j ) numerical 

and Av 2 and AV3 (Av x < 0.0002 
throughout) are plotted in sketch (i). 

A second modification to 
Treanor r s method would be to 

(l) Use the same differencing 
scheme 



Sketch (g).- Treanor 1 s modified method (B13) 
with step-size control (B12) applied to 
equations (Bll). 


5 Compute either by equation (B2) or (B 6), the results being the same. 
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(2) Set all elements of P equal to 8/h at every step. 

Such a method also has a real stability boundary |AH| = 0.5* The results of 
applying it to equation (Bll) are shown in sketches (j)> (k ) } and ( Z) . When 
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Sketch (h) . - Values of P calculated using 
Treanor* s modified method (B13) * 



Sketch (i).- Errors in v 2 and v 3 using Treanor T s 
modified method (B13) . 



Sketch (j).- Treanor 1 s modified method (Bl4) 
with step -size control (B12) applied to 
equations (Bll) . 



Sketch (k) . - Value of P 3 calculated using 

Treanor* s modified method (B14) . Variations 
of Pi and P 2 are nearly identical with P 3 . 





1 l I l 

,6 .8 t.O 1.2 

S 


Sketch (l).- Error in v 3 using Treanor 1 s 
modified method (Blij-) . 
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combined with the step-size control (B12), and started with As = 0.005, this 
method at once chose a step size equal to 0.01 for which it is marginally 
stable. The eigenvalue -500 was detected and -1000 excluded as long as this 
step size was used, see sketch (k) . Eventually, however, the step was doubled 
to the unstable size 0.01. In a few steps the eigenvalue -1000 emerged (as 
indicated by sketch (k) ) , and the error phenomenon typical of this kind of 
step control ensued. 

Certainly the curves for w in sketches (d), (g), and (j) could be made 
much smoother with more refined techniques for modifying the step size. How- 
ever, two things must be borne in mind: First, such refinements reduce the 

total distance the methods advance a solution in a given number of steps; 
second, the implicit method requires no such refinements to improve the 
smoothness . 
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APPENDIX C 


STABILITY CRITERION FOR TREANOR'S METHOD APPLIED TO v' = Av 


Treanor's method can be written in predictor -corrector form as 


(1) 1 , 

V(i/2) = V n + 2 hV ’ 


n 


(Cla) 


v 


( 2 ) 

n+( 1 / 2 ) 


= v n + 77 hv 


I vJ 1 ) 


2 uv n+(i/2) 


(Clb) 


( 3 ) 

Vl ' + h 


2 v i+(U ) F2 + T i'u / 2 ) P1,F2 + T -' (Fl ' ***> 


(Clc) 


v n+i = v n + hv n' F i + h8 3( Pv n + v n ') + hb £ 


C 1 ) (i)' 

^+( 1 / 2 ) + V n+(i/2) 


+ h5 2 


Py 


( 2 ) . ,.( 2 ) 


n+ ( 1 / 2 ) + v n +(i/2) 


+ h6i 


„ ( 3 ) ( 3 )' 

Pv^ 1 + Y v 

n+i n+i 


(Cld) 


where 


-Ph n 

i? = e ~ 1 

Fl -Ph ’ 


F _ e~ Ph - 1 + Ph 

F 2 - , . o > 


(Ph)' 


F 3 = 


e“ Ph - 1 + Ph - | (Ph) 2 


-(Ph)' 


(C2) 


and 

= -F2 + 4F3 

6 2 = 2(F 2 - 2F 3 ) \ (C3) 

5 3 = 4F 3 - 3 P 2 

The stability of equations (Cl) is determined by applying them to v' = Av. 

The resulting matrix equation 
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E 1/2 

0 

0 

- ( 1+ 1 *“) 

r 

- | AhE l/2 

E 1 / 2 

0 

-1 


-APh 2 F 2 E 1/2 

-2AhF 2 E 1/2 

E 

-l-Ah(Fi-2F 2 ) 

1 

V ] 

_2(P+A)h5 2 E l/2 

-(P+A)h& 2 E l/2 

-(P+A)hSiE 

E-(l+AhFi) - ( P+A ) h& 3 

Jf] 


'n 

n 

( 3 ) 


= 0 


(C4) 


where E = exp[h(d/dx)] has a characteristic equation with only one nonzero 
root. It is given by the equation 


E = 1 + AhFi + A(A + P)h 2 [2F 2 - 4F 3 + (4F 3 - F 2 )(Fi + PhF 2 )] 

+ A 2 (A + P)h 3 [2F 2 - 4F 3 + 2F 2 (4F 3 - F 2 )(2 + Ph) ] 

+ | A 3 (A + P)h 4 F 2 ( 4F 3 - F 2 ) (C5) 


The combinations of real positive Ph and negative Ah which make the right- 
hand side of equation (C5) equal to unity form the stability boundaries shown 
in figure 1. One can easily show 


lim PhFi = 1 , PbF 2 = 1 , PhF 3 = 1/2 

Ph-+oo 

Thus, for large values of Ph, equation (C5) reduces to 


E = 1 + Ah + | A 2 h 2 

which is the characteristic equation for the second -order, Runge-Kutta method. 
This gives the asymptotic value of the left stability boundary shown in 
figure 1. 
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